---
title: "EDA of NPD.no "
author: "Kristijan Bakaric"
date: "2018-10-21"
output:   
  flexdashboard::flex_dashboard:
    theme: journal
    social: menu
    source_code: embed
---


```{r setup, include=FALSE}
# prep workspace
library(data.table)# convenient reading function
library(tidyverse) # all you will ever need (almost all)
library(leaflet)   # interative mapping
library(DT)        # interactive tables
library(crosstalk) # inter-widget interactivity
library(sf)        # spatial data
library(mapview)   # spatial data eda
library(janitor)   # cleaning data
library(plotly)    # interactive plots
library(readxl)    # read excel files
```

Intro {data-icon="fa-picture-o"}
===================================== 

```{r, out.width = "100%"}
library(knitr)
library(png)
include_graphics("https://wallpapercave.com/wp/wp2681827.jpg")

# 

```


E&P of HC {data-icon="fa-info-circle"}
===================================== 

Sidebar {.sidebar}
----------------------------------------------------------------------

```{r}
# load summary wellbore table
wellbore_count_NCS <- read_excel("data/tabular/wellbore_count_figure_entry.xls", 
    sheet = "WELLS_KBAKA")
# gather the data for plotting
wellbore_count_NCS_gathered <- wellbore_count_NCS %>% 
  gather(key = YEAR, value = VALUE, -AREA, -TYPE, -PURPOSE)

```

```{r}
# total count of the wellbores in  Norway
wells_count <- wellbore_count_NCS_gathered %>% 
  summarise(count = sum(VALUE))

```




Exploration and production of oil and gas in Norway has a long history, starting
in **`r wellbore_count_NCS_gathered %>% summarize(min(YEAR))`**.

In total there has been drilled **`r wells_count`** wells, both for exploration and 
development purposes.

All wildcat and production wells drilled on the Norwegian shelf since 1966 have been registered in the Fact Pages with wellbore designation and well ID. Scientific wells and shallow drilling activity are not included in this list.

"The companies should not compete for raw data. Everyone should have access to basic geological data and rather compete to develop the best interpretation," says Robert Williams. He is a geologist and coordinator for the Core Store.

In total, 13,623,317 metres of production wells and 5,065,273 metres of exploration wells have been drilled.


Column {vertical_layout: scroll}
-----------------------------------------------------------------------




```{r well count plot per AREA}
# create a plot "Numbers of Wells Drilled in Norway by Geography"
wells_per_region <- ggplot(wellbore_count_NCS_gathered, aes(YEAR, VALUE, fill = AREA)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, 
                                   hjust = 1), 
                                   legend.position="bottom") +
  labs(x = "Year Drilled", 
       y = "Well Count",
       title = "Numbers of Wells Drilled in Norway by Geography",
       caption = "source: www.npd.no") 

ggplotly(wells_per_region, height = 600, width = 1200)
```





```{r well count plot per TYPE}
# create a plot "Numbers of Wells Drilled in Norway by Well Type"
wells_type <- ggplot(wellbore_count_NCS_gathered, aes(YEAR, VALUE, fill = TYPE)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position="bottom") +
  labs(x ="Year Drilled", y="Well Count",
       title ="Numbers of Wells Drilled in Norway by Well Type",
       caption ="source: www.npd.no") 

ggplotly(wells_type, height = 600, width = 1200)

```





```{r well count plot per PURPOSE}
# create a plot "Numbers of Wells Drilled in Norway by Well Purpose"
wells_purpose <- ggplot(wellbore_count_NCS_gathered, aes(YEAR, VALUE, fill = PURPOSE)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = "bottom") +
  labs(x = "Year Drilled", y = "Well Count",
       title = "Numbers of Wells Drilled in Norway by Well Purpose",
       caption = "source: www.npd.no")  

ggplotly(wells_purpose, height = 600, width=1200)
```


```{r load well data, include=FALSE}

# load ESRI's spatial data (point data)
wells_shp <- st_read("./data/spatial/v_geo_wlbpoint.shp") %>% 
  select(WBNAME, FIELD, geometry) %>% 
  mutate(centroid = st_point_on_surface(geometry)) # force centroid within a polygon

wells_shp_coords <- wells_shp$centroid %>%
  # extract the coordinates of these points as a matrix
  st_coordinates()
# insert centroid long and lat fields as attributes of polygons
wells_shp$lng <- wells_shp_coords[ , 1]
wells_shp$lat <- wells_shp_coords[ , 2]


wells_lithostratigraphy <- fread("http://factpages.npd.no/ReportServer?/FactPages/TableView/wellbore_formation_top&rs:Command=Render&rc:Toolbar=false&rc:Parameters=f&rs:Format=CSV&Top100=false&IpAddress=23.226.133.134&CultureCode=en", encoding = "UTF-8") %>% 
  filter(lsuLevel == "FORMATION") %>% 
  select(wlbName, lsuName, lsuTopDepth, lsuBottomDepth, lsuNameParent) 
  

wells_enriched <- wells_shp %>% 
  right_join(wells_lithostratigraphy, by = c("WBNAME" = "wlbName"))

length(unique(wells_enriched$WBNAME))

sd_wells <- SharedData$new(wells_enriched)


```


Wells - lithostrat {data-icon="fa-tint"}
=====================================

Well Filters {.sidebar}
-------------------------------------

```{r}
filter_checkbox(id = "lsuName" ,
                label = "Formation Name",
                sharedData = sd_wells,
                group = ~lsuName)

# setDefaultSelection(sd$groupName(), row.names(mtcars %>% filter(mpg > 30)
```

Row {.tabset .tabset-fade}
------------------------------------- 

### Map 

```{r map box}

# leaflet(sd_wells) %>% 
#   addTiles(providers$CartoDB.Positron) %>% 
#   addCircles(~lng, ~lat)
#   # label(~WBNAME)



sd_wells %>%
  plot_mapbox(lat = ~lat, lon = ~lng,text = ~WBNAME, color = ~lsuName,
              mode = 'scattermapbox', hoverinfo='text') %>%
  layout(title = 'Lithostratigraphy of the Wells',
         font = list(color='white'),
         plot_bgcolor = '#191A1A', paper_bgcolor = '#191A1A',
         mapbox = list(style = 'dark'),
         legend = list(orientation = 'h',
                       font = list(size = 8)),
         margin = list(l = 50, r = 50,
                       b = 100, t = 100,
                       pad = 4),
         autosize = F, width = 1400, height = 800, margin = margin)

```


### Table with Litostratigraphy

```{r data table}

DT::datatable(sd_wells,
  filter = "top",  # allows filtering on each column
  extensions = c(
    "Buttons"  # add download buttons, etc
    # "Scroller"  # for scrolling down the rows rather than pagination
  ),
  rownames = FALSE,  # remove rownames
  style = "bootstrap",
  class = "compact",
  width = "100%",
  options = list(
    dom = "Blrtip",  # specify content (search box, etc)
    deferRender = TRUE,
    scrollY = 300,
    scroller = TRUE,
    columnDefs = list(
      list(
        visible = FALSE
         
      )
    ),
    buttons = list(
      I("colvis"),  # turn columns on and off
      "csv",  # download as .csv
      "excel"  # download as .xlsx
    )
  )
)

  # colnames = c(
  #   "Discovery Name" = "DISCNAME",
  #   "Field Name" = "FIELDNAME",
  #   "Included in Field" = "INC_IN_FLD",
  #   "Discovery Status" = "DSCACTSTAT",
  #   "Company Long Name" = "CMP_LONG_NAME",
  #   "Discovery Status" = "DSC_CURRENT_ACTIVITY_STATUS",
  #   "HC type" = "MAP_DSC_HC",
  #   "Wellbore Name" = "WLB_NAME",
  #   "Discovery Year" = "DSC_DISCOVERY_YEAR",
  #   "Discovery Owner Type" = "DSC_OWNER_KIND",
  #   "Discovery Owner Name" = "DSC_OWNER_NAME"
  # )


```


```{r load discoveries data, include=FALSE}
# load the table containing overview of discoveries from npd.no and clean the names
discoveries_overview <- fread("http://factpages.npd.no/ReportServer?/FactPages/TableView/discovery&rs:Command=Render&rc:Toolbar=false&rc:Parameters=f&rs:Format=CSV&Top100=false&IpAddress=193.71.188.220&CultureCode=en", encoding = "UTF-8") %>% 
  clean_names(case = "screaming_snake")

# load polygons representing discoveries (local shp files) :(
discoveries_spatial <- st_read("./data/spatial/v_geo_dscarea.shp") %>% 
  mutate(centroid = st_point_on_surface(geometry)) # force a central point in polygon

# extract coordinates from a centroid points since I could not get the polygons 
# to work with the crosstalk package (i.e. with brush filtering)
discoveries_spatial_coords <- discoveries_spatial$centroid %>%
  # extract the coordinates of these points as a matrix
  st_coordinates()
# insert centroid long and lat fields as attributes of polygons
discoveries_spatial$long <- discoveries_spatial_coords[ , 1]
discoveries_spatial$lat <- discoveries_spatial_coords[ , 2]

# enrich shp discoveries with more attributes
discoveries_joined <- discoveries_spatial %>% 
  left_join(discoveries_overview, by = c("IDDISCOVER" = "DSC_NPDID_DISCOVERY")) %>% 
  select(DISCNAME, FIELDNAME, DSCACTSTAT, INC_IN_FLD, MAP_DSC_HC, 
         CMP_LONG_NAME, DSC_CURRENT_ACTIVITY_STATUS, WLB_NAME,
         NMA_NAME, DSC_DISCOVERY_YEAR, DSC_DATE_FROM_INCL_IN_FIELD,  DSC_OWNER_KIND,  DSC_OWNER_NAME, LABEL, geometry, centroid, lat, long) 

# create a shared object for interactivity (crosstalk package)
sd <- SharedData$new(discoveries_joined)

# workaround if ever needed to use sp package
# discoveries_joined <- discoveries_joined %>% as('Spatial')
# sd_map <- SharedData$new(discoveries_joined)
# sd_df <- SharedData$new(as.data.frame(discoveries_joined@data), group = sd_map$groupName())

```




Production {data-icon="fa-line-chart"}
===================================== 

Column Inputs {.sidebar}
-----------------------------------------------------------------------

**Cumulative production** is defined as the gross amount of oil and gas production from an oil reservoir over a particular period of time of the life of a well. In general, cumulative production is an oil and gas industry term related to an oil well, a basin or an oil field.
In this case you can inspect cummulative curve of an hydrocarbon field by clicking on a 
field polygon

In the **Interactive map** you can:

* click to grab and drag the map around
* zoom with the '+' and '-' buttons or with your mouse's scroll wheel
* click a field polygon to reveal a popup with the cummulative production chart



Column {vertical_layout: scroll}
-----------------------------------------------------------------------

```{r yearly production processing, message=FALSE, warning=FALSE, echo=FALSE, include = FALSE}

# load data about yearly hydrocarbon production in Norway

fields_yearly_production <- fread("http://hotell.difi.no/download/npd/field/production-yearly-by-field?download", encoding = "UTF-8") %>% 
    mutate_at(vars(ends_with("Sm3")), funs(as.numeric)) 

# 
fields_yearly_production_OeNetMillSm3 <- fields_yearly_production %>% 
  select(prfInformationCarrier, prfNpdidInformationCarrier, prfYear, prfPrdOilNetMillSm3,
         prfPrdGasNetBillSm3, prfPrdNGLNetMillSm3, 
         prfPrdCondensateNetMillSm3, prfPrdProducedWaterInFieldMillSm3) %>% 
  mutate(prfPrdGasNetMillSm3 = (prfPrdGasNetBillSm3 )) %>% 
  select(-prfPrdGasNetBillSm3) %>% 
  gather(key = HC_PHASE_MillSm3, value = value, -c(prfInformationCarrier, prfNpdidInformationCarrier, prfYear)) %>% 
  arrange(prfInformationCarrier, prfYear, HC_PHASE_MillSm3) %>% 
  group_by(prfInformationCarrier, prfNpdidInformationCarrier, HC_PHASE_MillSm3) %>%
  mutate(cumsumOEprod = cumsum(value)) %>% 
  group_by(prfInformationCarrier, prfNpdidInformationCarrier) %>% 
  nest() # nest the data per hydrocarbon field in tidy format



# fields_yearly_production_OeNetMillSm3[42,]$data

# generate plots for yearly production per field in a tidy format (as a plot column)
fields_yearly_production_OeNetMillSm3 <- fields_yearly_production_OeNetMillSm3 %>%
  mutate(plot = map2(data, prfInformationCarrier,  ~ggplot(data = .x) +
       geom_col(aes(x = prfYear, y = cumsumOEprod, fill = HC_PHASE_MillSm3 )) +
       ggtitle(.y) +
       theme_bw() +
       theme(legend.position="bottom") +
       ylab("MillSm3") +
       xlab("Year") +
       scale_fill_manual(values = alpha(c("pink", "red","black","green","blue"), .4)) +
       theme(axis.text.x = element_text(angle = 90, hjust = 1))))



# fields_yearly_production_OeNetMillSm3[22:50, ]$plot

# load polygon field outlines (shp files)
fields_shp <- st_read("./data/spatial/v_geo_fldarea.shp") %>% 
    select(IDFIELD, FIELDNAME, IDDISCOVER, DISCNAME, geometry) %>% 
    filter(!is.na(FIELDNAME)) %>% 
    mutate(centroid = st_point_on_surface(geometry))
  

# enrich the field data with more attributes
fields_shp_enriched <- fields_shp %>% 
  left_join(fields_yearly_production_OeNetMillSm3, by = c("IDFIELD" = "prfNpdidInformationCarrier")) %>% 
  filter(!is.na(prfInformationCarrier))

```



```{r yearly production map, mapview}
# explore the yearly field production per year in a combination of a map + pop-up chart

mapviewOptions(leafletWidth = "1600px",leafletHeight = "1000px")

mapview(fields_shp_enriched, 
        popup = popupGraph(fields_shp_enriched$plot, 
                           height = 600, 
                           width = 800),
        label = fields_shp_enriched$DISCNAME)

```



Discoveries {data-icon="ion-stats-bars"}
=====================================  

Column Inputs {.sidebar}
-------------------------------------

### Filters

```{r filters}
filter_select(
  id = "DISCNAME",
  label = "Discovery Name",
  sharedData = sd,
  group = ~DISCNAME
)

  filter_checkbox(
    id = "MAP_DSC_HC",
    label = "Discovery Type",
    sharedData = sd,
    group = ~MAP_DSC_HC
  )

filter_checkbox(
    id = "NMA_NAME",
    label = "Region Name",
    sharedData = sd,
    group = ~NMA_NAME
  )

filter_checkbox(
    id = "DSCACTSTAT",
    label = "Discovery Status",
    sharedData = sd,
    group = ~DSCACTSTAT
  )

filter_slider(
    id = "DSC_DISCOVERY_YEAR",
    label = "Discovery Year",
    sharedData = sd,
    column = ~DSC_DISCOVERY_YEAR,
    step = 1,
    round = TRUE,
    sep = "",
    ticks = FALSE
  )

  

```

Row {data-height=550}
-------------------------------------

### Datatable


    
```{r datatable}
sd %>% 
  DT::datatable(
    filter = "top",  # allows filtering on each column
    extensions = c(
      "Buttons"  # add download buttons, etc
      # "Scroller"  # for scrolling down the rows rather than pagination
    ),
    rownames = FALSE,  # remove rownames
    style = "bootstrap",
    class = "compact",
    width = "100%",
    options = list(
      dom = "Blrtip",  # specify content (search box, etc)
      deferRender = TRUE,
      scrollY = 300,
      scroller = TRUE,
      columnDefs = list(
        list(
          visible = FALSE,
          targets = c(14:17)
        )
      ), 
      buttons = list(
        I("colvis"),  # turn columns on and off
        "csv",  # download as .csv
        "excel"  # download as .xlsx
      )
    ),
    colnames = c(
      "Discovery Name" = "DISCNAME",
      "Field Name" = "FIELDNAME",
      "Included in Field" = "INC_IN_FLD",
      "Discovery Status" = "DSCACTSTAT",
      "Company Long Name" = "CMP_LONG_NAME",
      "Discovery Status" = "DSC_CURRENT_ACTIVITY_STATUS",
      "HC type" = "MAP_DSC_HC",
      "Wellbore Name" = "WLB_NAME",
      "Discovery Year" = "DSC_DISCOVERY_YEAR",
      "Discovery Owner Type" = "DSC_OWNER_KIND",
      "Discovery Owner Name" = "DSC_OWNER_NAME"
    )
  )
```

Row {data-height=550}
-------------------------------------
    
### Interactive map
    
```{r map}


getColor <- case_when(discoveries_joined$MAP_DSC_HC == "OIL/GAS" ~ "orange",
          discoveries_joined$MAP_DSC_HC == "GAS" ~ "red",
          discoveries_joined$MAP_DSC_HC == "GAS/CONDENSATE" ~ "pink",
          discoveries_joined$MAP_DSC_HC == "OIL" ~ "black")

icons <- awesomeIcons(
  icon = 'flame',
  iconColor = 'yellow',
  library = 'ion',
  markerColor = getColor
)

# set your own colors manually:
pal <- colorFactor(
  palette = c('red', 'pink', 'green'),
  domain = discoveries_joined$MAP_DSC_HC
)

leaflet(sd) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
  label =  ~LABEL,
  fillColor = ~pal(MAP_DSC_HC),
  weight = 2,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE))  %>%
  addAwesomeMarkers(~long, 
             ~lat,
             popup = ~as.character(LABEL),
       icon = icons
    ) %>% 
  addMiniMap(
    tiles = providers$CartoDB.Positron,
    toggleDisplay = TRUE)


# leaflet(df.20) %>% addTiles() %>%
#   addAwesomeMarkers(~long, ~lat, icon=icons, label=~as.character(mag))

```



Investment vs Production {data-orientation=rows data-icon="fa-money"}
===================================== 


```{r}
ProdFieldYearSellable <- fread("http://hotell.difi.no/download/npd/field/production-yearly-by-field", encoding = "UTF-8")
ProdFieldYearSellable <- ProdFieldYearSellable %>% 
  mutate_at(vars(ends_with("Sm3")), funs(as.numeric))

InvestmentsFieldYear <- fread("http://hotell.difi.no/download/npd/investments/yearly-by-field", encoding = "UTF-8")

FieldsDescriptions <- fread("http://hotell.difi.no/download/npd/field/overview", encoding = "UTF-8")

OperatorsByField <-  fread("http://hotell.difi.no/download/npd/field/operators", encoding = "UTF-8")

```

```{r}

InvestmentsFieldYear <- InvestmentsFieldYear %>% 
  group_by(prfInformationCarrier) %>% 
  mutate(cumulativeInvestmentsMillNOK = cumsum(prfInvestmentsMillNOK))

ProdFieldYearSellable <- ProdFieldYearSellable %>% 
  group_by(prfInformationCarrier) %>% 
  mutate(cumulativeprfPrdOeNetMillSm3 = cumsum(prfPrdOeNetMillSm3))

ProdPlusInvestmentYear <- InvestmentsFieldYear %>% 
  full_join(ProdFieldYearSellable, by = c("prfInformationCarrier", "prfYear"))  %>% 
  left_join(subset(FieldsDescriptions, select = c("fldName","fldMainArea")), by = c("prfInformationCarrier" = "fldName")) 

 ProdPlusInvestmentYear <- ProdPlusInvestmentYear %>% 
   na.omit(prfInvestmentsMillNOK) %>% 
   na.omit(prfPrdOeNetMillSm3)

```

### North Sea

```{r}
library(plotly)

NorthSea <- ggplot(subset(ProdPlusInvestmentYear, fldMainArea == "North sea"), aes(x = cumulativeInvestmentsMillNOK, y = cumulativeprfPrdOeNetMillSm3, color = fldMainArea))+
  geom_point(aes(size = cumulativeInvestmentsMillNOK, frame = prfYear, ids = prfInformationCarrier))
NorthSea <- ggplotly(NorthSea)
NorthSea <- NorthSea %>% 
  animation_opts(
    3000, easing = "elastic", redraw = FALSE
  )

NorwegianSea <- ggplot(subset(ProdPlusInvestmentYear, fldMainArea == "Norwegian sea"), aes(x = cumulativeInvestmentsMillNOK, y = cumulativeprfPrdOeNetMillSm3, color = fldMainArea)) + 
  geom_point(aes(size = cumulativeInvestmentsMillNOK, frame = prfYear, ids = prfInformationCarrier))
NorwegianSea <- ggplotly(NorwegianSea)
NorwegianSea <- NorwegianSea %>% 
  animation_opts(
    3000, easing = "elastic", redraw = FALSE
  )

BarentsSea <- ggplot(subset(ProdPlusInvestmentYear, fldMainArea == "Barents sea"), aes(x = cumulativeInvestmentsMillNOK, y = cumulativeprfPrdOeNetMillSm3, color = fldMainArea))+
  geom_point(aes(frame = prfYear, ids = prfInformationCarrier))
BarentsSea <- ggplotly(BarentsSea)
BarentsSea <- BarentsSea %>% 
  animation_opts(
    3000, easing = "elastic", redraw = FALSE
  )

NorthSea
```

### Norwegian Sea

```{r}
NorwegianSea
```

### Barents Sea

```{r}
BarentsSea
```

Information {data-orientation=rows data-icon="fa-info-circle"}
===================================== 

R credits {.sidebar}
-------------------------------------

**Credits:**    
The following packages were used: [flexdashboard](https://cran.r-project.org/web/packages/flexdashboard/index.html), [rgdal](https://cran.r-project.org/web/packages/rgdal/index.html), [dplyr](https://cran.r-project.org/web/packages/dplyr/index.html), [leaflet](https://cran.r-project.org/web/packages/leaflet/index.html), [d3scatter](https://github.com/jcheng5/d3scatter), [DT](https://cran.r-project.org/web/packages/DT/index.html), [gapminder](https://cran.r-project.org/web/packages/gapminder/index.html), [countrycode](https://cran.r-project.org/web/packages/countrycode/index.html),
[rworldmap](https://cran.r-project.org/web/packages/rworldmap/index.html), and [spatialEco](https://cran.r-project.org/web/packages/spatialEco/index.html).


Data Sources {.sidebar}
-------------------------------------

text